---
title: "Data Science Transport – Second Assignment – Group 12"
output: 
  flexdashboard::flex_dashboard:
    orientation: columns
    vertical_layout: fill
    source_code: embed
---

```{r setup, include=FALSE}
library(flexdashboard)
library(gtfsrouter)
library(tidyverse)
library(tidytransit)
library(sf)
library(tmap)
library(units)
library(RColorBrewer)
tmap_mode("view")
```

last time {data-icon="fa-hourglass-half"}
=====================================
TEMPLATE BULLSHIT

TO EDIT

ADD DATA SOURCES!!!!!!!!!!!!!!!

- test1
- test2

### Chart 0
```{r}

library(leaflet)
leaflet() %>%
  addTiles() %>%
  addMarkers(lng=174.768, lat=-36.852, popup="The birthplace of R")
```

### Chart 1

```{r}
# 1. Plot the dots themselves
```

### Cars

```{r}
knitr::kable(mtcars)
```

isochrones {data-icon="fa-expand-arrows-alt"}
=====================================

```{r, include = FALSE}
##############################################################
#
#   READ GTFS DATA
#
##############################################################
# set work directions
setwd_gtfs <- function(){setwd("~/Documents/Uni/Master/DataScienceTransport/data/vbb-gtfs")}
setwd_data <- function(){setwd("~/Documents/Uni/Master/DataScienceTransport/data")}
setwd_work <- function(){setwd("~/Documents/Uni/Master/DataScienceTransport/assignment_2")}

setwd_work

# read gtfs data for monday
file <- file.path("~/Documents/Uni/Master/DataScienceTransport/data/vbb-gtfs/2020-12_2020-12-28.zip")
gtfs <- extract_gtfs(file) %>% gtfs_timetable(day = 2)

##############################################################
#
#   SET TIMES
#
##############################################################
start_time <- 7 * 3600 + 1200
end_time <- 8 * 3600

# create isochrone
# ic <- gtfs_isochrone (gtfs,
#                       from = from,
#                       start_time = start_time,
#                       end_time = end_time)

##############################################################
#
#   CREATE STOPS SF OBJECT
#
##############################################################
stops <- st_as_sf(gtfs$stops,
                   coords = c("stop_lon", "stop_lat"),
                   crs = 4326) %>% 
  st_transform(25833)

##############################################################
#
#   SHAPE DISTRICTS NEW (+ area)
#
##############################################################
setwd_data()
shape_districts_new <- read_sf(dsn = "LOR_SHP_2019-1", layer = "Planungsraum_EPSG_25833")
setwd_work()

shape_districts_new <- shape_districts_new %>% 
  group_by(BEZIRK) %>% 
  summarise() %>% 
  filter(!is.na(BEZIRK)) %>% 
  rename(NAME = BEZIRK) %>% 
  mutate(AREA = st_area(geometry)) %>% 
  select(NAME, AREA, everything()) %>% 
  mutate(AREA = (AREA / 1000000) * as_units("km2"))

# setting crs of polygons
st_crs(shape_districts_new$geometry) <- 25833

shape_berlin <- st_union(shape_districts_new)

##############################################################
#
#   SPECIFIC SHAPES AND STOPS
#
##############################################################

stops_in_berlin <- stops %>% 
  mutate(inside_berlin = st_within( geometry, shape_berlin )) %>% 
  mutate(inside_berlin = !is.na( as.numeric( inside_berlin ))) %>% 
  filter(inside_berlin == TRUE)

# get isochrone area
# ic = gtfs_isochrone (gtfs,
#                      from = "Berlin, Sowjetisches Ehrenmal",
#                      #from_is_id = TRUE,
#                      start_time = start_time,
#                      end_time = end_time)$hull$area

##############################################################
#
#   CALCULATE ISOCHRONES
#
##############################################################
# # the following code calculates the isochrones (inefficent, ~ 15h)
# # instead of running the code, we read in the pre-calculated file

# stops_ic_area <- vector(mode = "double")
# 
# # create isochrone areas for stops in 50 minutes
# for (stop_name in stops$stop_name){
# 
#   tryCatch( {
#     ic_area <- gtfs_isochrone (gtfs,
#                                from = stop_name,
#                                #from_is_id = TRUE,
#                                start_time = start_time,
#                                end_time = end_time)$hull$area
#     if(is.null(ic_area)) {
#       stops_ic_area <<- rbind(stops_ic_area, 0)
#       print(paste(stop_name, ": ", ic_area, "!!!!!!!!!!"))
#     } else {
#       stops_ic_area <<- rbind(stops_ic_area, ic_area)
#       print(paste(stop_name, ": ", ic_area))
#     }
#     },
#     error = function(e) {
#       stops_ic_area <<- rbind(stops_ic_area, 0)
#       print(paste("ERROR!!!", stop_name))
#       }
#     )
# }
# 
# ##############################################################
# #
# #   CLEANING
# #
# ##############################################################
# 
# # merge and clean
# # https://r-spatial.github.io/sf/reference/bind.html
# # https://cran.r-project.org/web/packages/units/vignettes/units.html
# rownames(stops_ic_area) <- NULL
# stops_area <- st_sf(data.frame(stops, stops_ic_area / 1000000)) %>%
#   rename(ic_area = stops_ic_area.1e.06,
#          id = stop_id,
#          name = stop_name,
#          parent = parent_station) %>% 
#   select(id, name, parent, ic_area) %>% 
#   mutate(ic_area = ic_area * as_units("km2"))
# 
# # save
# # https://r-spatial.github.io/sf/reference/st_write.html
# st_write(stops_area, "output_stops_ic_area.shp")
stops_area <- st_read("output_stops_ic_area.shp")

# more cleaning for plot
# https://dplyr.tidyverse.org/reference/distinct.html
stops_area = 
  stops_area %>% 
  select(name, ic_area) %>% 
  distinct(name, .keep_all = TRUE)

stops_area_berlin <- stops_area %>% 
  mutate(inside_berlin = st_within( geometry, shape_berlin )) %>% 
  mutate(inside_berlin = !is.na( as.numeric( inside_berlin ))) %>% 
  filter(inside_berlin == TRUE) %>% 
  select(-inside_berlin) %>% 
  mutate(id = paste(name, ": ", round(ic_area)))
```

Column {data-width=100}
-------------------------------------
    
### Approach {data-height=200}

The [gtfsrouter](https://atfutures.github.io/gtfs-router/) package allows us to calculate all stations reachable within a specified time period from a nominated station ([isochrones](https://atfutures.github.io/gtfs-router/reference/gtfs_isochrone.html)). We use the hull polygon as an indicator for the city area reachable.

See below for a isochrone plot of a 30min Monday work–home trip starting from Helmholtzstraße at 18:00.

We try to simulate a home–work trip in the morning rush hour arriving at 08:00. Assuming a five minute walk at beginning and end of the trip, we limit the GTFS data to Monday (2021-01-18) and calculate the isochrones for a 40min time range. By that we want to get all reachable stations by local transport (< 60min, [PBefG §8 (1)](https://www.gesetze-im-internet.de/pbefg/__8.html)) and rate the station inside the network.

- reachablity of possible work places
- calculation with all transport agencies
-- domination of the hull area in 40 min time range by reachable fanning out long-distance transport
- no weighting of the reachable area

### 30min IC home trip (18:00) from Helmholtzstr. {data-height=100}
    
```{r}
ic_einstein <- gtfs_isochrone(gtfs,
                              from = "Berlin, Helmholtzstr.",
                              start_time = 18 * 3600,
                              end_time = 18 * 3600 + 1800)

tm_basemap(leaflet::providers$OpenStreetMap.DE) +
  tm_shape(ic_einstein$hull) + 
  tm_polygons(col = "red",
              alpha = 0.2,
              border.col = "red") +
  tm_shape(ic_einstein$routes) +
  tm_lines() +
  tm_shape(ic_einstein$end_points) +
  tm_dots(col = "red") + 
  tm_shape(ic_einstein$start_point) + 
  tm_dots(col = "green")
```

Column {data-width=300}
-------------------------------------
   
### area size of the hull enclosing the routed points

```{r}
##############################################################
#
#   PLOT
#
##############################################################

tm_shape(shape_districts_new) +
  tm_polygons(alpha = 0,
              popup.vars = c("area" = "AREA")) +
  tm_shape(stops_area_berlin) +
  tm_dots(col = "ic_area",
          id = "name",
          popup.vars = c("area" = "ic_area"),
          size = 0.07,
          border.lwd = 0.3,
          legend.hist = TRUE,
          n = 15,
          title = "isochrone area [km^2]") +
  tm_view(bbox = shape_berlin)
```


traveltimes to center {data-icon="fa-stopwatch"}
=====================================

```{r, include = FALSE}
##############################################################
#
#   SHAPE CENTER AREAS
#
##############################################################
# "Zentrentragender Stadtraum mit höchster / hoher Urbanität"
# of Zentrumsbereichskernen
# see page 39: https://www.stadtentwicklung.berlin.de/planen/stadtentwicklungsplanung/download/zentren/2011-07-31_StEP_Zentren3.pdf
# or page 45 (less detailed): https://www.stadtentwicklung.berlin.de/planen/stadtentwicklungsplanung/download/zentren/StEP_Zentren_2030.pdf
# recreated with QGis

shape_center <- read_sf(dsn = "shape_center_areas", layer = "center_areas") %>% 
  mutate(name = c("east", "west")) %>% 
  select(name)

shape_center_east <- shape_center %>% filter(name == "east")
shape_center_west <- shape_center %>% filter(name == "west")

##############################################################
#
#   READ GTFS DATA
#
##############################################################
# now we work with tidytransit
# calculation of shortest tt from all station to specific ones is more convinent

setwd_gtfs()
gtfs <- read_gtfs("2020-12_2020-12-28.zip")
setwd_work()

# http://tidytransit.r-transit.org/reference/filter_stop_times.html
stop_times_filtered <- filter_stop_times(gtfs, "2021-01-18", "06:00:00", "07:55:00")

##############################################################
#
#   GET STOPS
#
##############################################################
stops <- st_as_sf(gtfs$stops, coords = c("stop_lon", "stop_lat"), crs = 4326) %>%
  st_transform(25833) %>% 
  select(stop_name) %>%
  rename(name = stop_name) %>%
  distinct(name)

stops_berlin <- stops %>% 
  mutate(inside_berlin = st_within( geometry, shape_berlin )) %>% 
  mutate(inside_berlin = !is.na( as.numeric( inside_berlin ))) %>% 
  filter(inside_berlin == TRUE) %>% 
  select(name)

stops_center <- stops %>% 
  mutate(inside_center = st_within( geometry, shape_center )) %>% 
  mutate(inside_center = !is.na( as.numeric( inside_center ))) %>% 
  filter(inside_center == TRUE) %>% 
  select(name)

stops_center_east <- stops %>% 
  mutate(inside_center_east = st_within( geometry, shape_center_east )) %>% 
  mutate(inside_center_east = !is.na( as.numeric( inside_center_east ))) %>% 
  filter(inside_center_east == TRUE) %>% 
  select(name)

stops_center_west <- stops %>% 
  mutate(inside_center_west = st_within( geometry, shape_center_west )) %>% 
  mutate(inside_center_west = !is.na( as.numeric( inside_center_west ))) %>% 
  filter(inside_center_west == TRUE) %>% 
  select(name)

##############################################################
#
#   TT calculation
#
##############################################################
# what are the tt to the center areas?
# according to Nahverkehrsplan Berlin 2019-2023: ANlage 1 - Monitoringbericht (p. 12)
# standard: tt_max = 3600, n_transfer_max = 2, n_realise_stations = 0.95

tt <- travel_times(
  stop_times_filtered,
  stops_center$name,
  time_range = 5400,
  arrival = TRUE,
  max_transfers = 2,
  # max_departure_time = NULL,
  return_coords = TRUE,
  return_DT = FALSE
)

# clean it for plot
tt <- tt %>% 
  rename(from = from_stop_name,
         to = to_stop_name,
         tt = travel_time,
         departure = journey_departure_time,
         arrival = journey_arrival_time
         ) %>% 
  select(-c(from_stop_id, to_stop_id, to_stop_lat, to_stop_lon)) %>% 
  st_as_sf(coords = c("from_stop_lon", "from_stop_lat"),
           crs = 4326) %>% 
  st_transform(25833) %>% 
  mutate(tt = set_units(round(tt/60, 2), "min"))
```

Column {data-width=100}
-------------------------------------
    
### approach {data-height=400}

The [local transport plan](https://www.berlin.de/sen/uvk/verkehr/verkehrsplanung/oeffentlicher-personennahverkehr/nahverkehrsplan/) (p. 106) sets targets for the connectivity standards. Different categories of center areas (see [StEP](https://stadtentwicklung.berlin.de/planen/stadtentwicklungsplanung/de/zentren/zentren2030/index.shtml), p. 45) have to be reachable within a certain time and with a maximum number of transfer. This should hold for 95% of the stations.

Based on the GTFS data, we try to recreate the result of the monitoring ([NVP Anlage 1](https://www.berlin.de/sen/uvk/_assets/verkehr/verkehrsplanung/oeffentlicher-personennahverkehr/nahverkehrsplan/broschure_nvp_2019_anlage_1.pdf), p. 12), mentioning a degree of fulfillment of 99,7% for the central areas.

* destination:
  + City West (Zoo/ Kurfürstendamm)
  + Mitte (Potsdamer Platz/ Alexanderplatz)
* max. tt: 60min
* max. transfers: 2

From a more detailed [display of the area](https://www.stadtentwicklung.berlin.de/planen/stadtentwicklungsplanung/download/zentren/2011-07-31_StEP_Zentren3.pdf) (p. 39), we create a shape file enclosing the associated stations. The [tidytransit](https://tidytransit.r-transit.org/) package let's us [calculate](https://tidytransit.r-transit.org/reference/travel_times.html) the shortest travel time for all stations to any of a specified set of stations. For that `arrival` has to be set to `TRUE`.

We try to simulate a home–work trip in the morning rush hour arriving at 08:00. Assuming five minute walking at beginning and end of the trip, we limit the GTFS data to Monday (2021-01-18) and calculate the shortest travel times in a 90min time range ending at 07:55.

### Percent of all berlin stations fullfill the connectivity standard according to the GTFS data {data-height=200}

```{r}
##############################################################
#
#   DEGREE OF FULLFILMENT
#
##############################################################

n_of_stations <- tt %>%
  mutate(inside_berlin = st_within( geometry, shape_berlin )) %>% 
  mutate(inside_berlin = !is.na( as.numeric( inside_berlin ))) %>% 
  filter(inside_berlin == TRUE) %>% 
  mutate(outside_center = st_within( geometry, shape_center )) %>% 
  mutate(outside_center = is.na( as.numeric( outside_center ))) %>% 
  filter(outside_center == TRUE) %>%
  nrow()

n_of_stations_valid <- tt %>% 
  mutate(inside_berlin = st_within( geometry, shape_berlin )) %>% 
  mutate(inside_berlin = !is.na( as.numeric( inside_berlin ))) %>% 
  filter(inside_berlin == TRUE) %>% 
  mutate(outside_center = st_within( geometry, shape_center )) %>% 
  mutate(outside_center = is.na( as.numeric( outside_center ))) %>% 
  filter(outside_center == TRUE) %>%
  filter(tt <= 60 * as_units("min")) %>% 
  filter(transfers <= 2) %>% 
  nrow()

percent_stations_valid <- n_of_stations_valid / n_of_stations * 100
percent_stations_valid <- round(percent_stations_valid, 2)

valueBox(paste(percent_stations_valid, "%"), icon = "fa-crosshairs")
```
    

Column {data-width=300}
-------------------------------------
   
### shortest travel time to one of the stations inside City West or Mitte

```{r}
##############################################################
#
#   PLOT
#
##############################################################

# https://campus.datacamp.com/courses/visualizing-geospatial-data-in-r/raster-data-and-color?ex=9
rdylgn <- rev(brewer.pal(7, "RdYlGn"))

# https://leaflet-extras.github.io/leaflet-providers/preview/
# https://tlorusso.github.io/geodata_workshop/tmap_package
# https://www.rdocumentation.org/packages/tmap/versions/3.0/topics/tm_basemap
# https://rdrr.io/cran/tmap/man/tm_view.html
# https://leafletjs.com/reference-1.3.4.html#map-methods-for-modifying-map-state

tm_basemap(leaflet::providers$CartoDB.DarkMatter) +
  tm_shape(shape_districts_new) + 
  tm_polygons(alpha = 0,
              lwd = 1.5,
              border.col = "white",
              popup.vars = c("area" = "AREA")
              ) +
  tm_shape(shape_center) +
  tm_polygons(alpha = 0.2,
              col = "red",
              border.col = "red"
              ) + 
  tm_shape(tt) +
  tm_dots(col = "tt",
          style = "fixed",
          breaks = c(0, 10, 20, 30, 40, 50, 60, 120),
          labels = c("0 – 10", "10 – 20", "20 – 30", "30 – 40", "40 – 50", "50 – 60", "> 60"), 
          id = "from",
          palette = rdylgn,
          title = "traveltime [min]",
          popup.vars = c("to" = "to", 
                         "traveltime" = "tt",
                         "departure at" = "departure",
                         "arrival at" = "arrival",
                         "number of transfers" = "transfers")
          ) +
  tm_view(bbox = shape_center)
```


***

https://rstudio.github.io/leaflet/

- Interactive panning/zooming

- Compose maps using arbitrary combinations of map tiles, markers, polygons, lines, popups, and GeoJSON.

- Create maps right from the R console or RStudio

- Embed maps in knitr/R Markdown documents and Shiny apps

- Easily render Spatial objects from the sp package, or data frames with latitude/longitude columns

- Use map bounds and mouse events to drive Shiny logic